Statistical analysis using causal inference

Agalychnis lemur

Author

Marcelo Araya-Salas, PhD

Purpose

  • Determine the adjustment sets that allow to infer a causal effect of environmental variables on vocal activity
  • Evaluate causal effect of environmental factors on vocal activity of A. lemur with bayesian regression models

Analysis flowchart

flowchart
  A[Define DAG] --> B(24 hours previous rain) 
  A --> C(48 hours previous rain)
  B --> D(Define adjustment sets for each predictor)
  C --> D
  D --> E(Run all models satisfying\nthe back door criterion)
  E --> F(Average posterior probabilities) 
  F --> G(Combine models in a single graph) 

style A fill:#44015466
style B fill:#3E4A894D
style C fill:#3E4A894D
style D fill:#26828E4D
style E fill:#6DCD594D
style F fill:#FDE7254D
style G fill:#31688E4D

 

Load packages

Code 
# Unload packages
out <- sapply(paste('package:', names(sessionInfo()$otherPkgs), sep = ""), function(x) try(detach(x, unload = FALSE, character.only = TRUE), silent = T))
  
pkgs <-
  c(
    "remotes",
    "viridis",
    "brms",
    "cowplot",
    "posterior",
    "readxl",
    "HDInterval",
    "kableExtra",
    "knitr",
    "ggdist",
    "ggplot2",
    "lunar",
    "cowplot",
    "maRce10/brmsish",
    "warbleR",
    "ohun",
    "dagitty",
    "ggdag",
    "tidybayes",
    "pbapply"
  )


# install/ load packages
sketchy::load_packages(packages = pkgs)
  
options("digits" = 6, "digits.secs" = 5, knitr.table.format = "html") 

# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE, tidy = TRUE)

# set working directory as project directory or one directory above,
opts_knit$set(root.dir = "..")

source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
source("~/Dropbox/R_package_testing/brmsish/R/check_rds_fits.R")
source("~/Dropbox/R_package_testing/brmsish/R/draw_extended_summary.R")

Custom functions

Code 
print_data_frame <- function(x){
    # print table as kable  
  kb <-kable(x, row.names = TRUE, digits = 3) 
    
  kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    
  print(kb)
}

adjustment_set_formulas <-
  function(dag,
           exposure,
           required_variable,
           outcome,
           effect = "total",
           type = "minimal",
           formula_parts = c(outcome),
           latent  = NULL,
           remove = NULL,
           plot = TRUE, 
           ...) {
    if (plot)
      gg <- ggdag_adjustment_set(
        .tdy_dag = tidy_dagitty(dag),
        exposure = exposure,
        outcome = outcome,
        ...
      ) + theme_dag()
    
    temp_set <-
      adjustmentSets(
        x = dag,
        exposure = exposure,
        outcome = outcome,
        effect = effect,
        type = type
      )
    
    
    form_set <- lapply(temp_set, function(x) {
      if (!is.null(remove))
        x <- x[!x %in% remove]
      form <-
        paste(
          formula_parts[1],
          " ~ ",
          exposure,
          " + ",
          paste(x, collapse =  " + "),
          if (length(formula_parts) == 2)
            formula_parts[2]
          else
            NULL
        )
      
      return(form)
    })
    
    form_set <- form_set[!duplicated(form_set)]
    
    if (!is.null(latent))
      for (i in latent)
        form_set <- form_set[!grepl(paste0(" ", i, " "), form_set)]
    
    # form_set <- sapply(form_set, as.formula)
    
    
    names(form_set) <- seq_along(form_set)
    
    # add formula as attribute
    attributes(form_set)$exposure.formula <- paste(formula_parts[1],
                                                   " ~ ",
                                                   exposure,
                                                   if (length(formula_parts) == 2)
                                                     formula_parts[2]
                                                   else
                                                     NULL)
    
    if (plot)
    return(gg) else
    return(form_set)
  }

# Define a function to remove special characters
remove_special_chars <- function(text) {
  # Replace special characters with hyphen
  cleaned_text <- gsub("[^a-zA-Z0-9]+", "-", text)
  # Remove leading and trailing hyphens
  cleaned_text <- gsub("^-+|-+$", "", cleaned_text)
  return(cleaned_text)
}

pa <- function(...)
  brms::posterior_average(...)

# to get average posterior values from models with different formulas
averaged_model <-
  function(formulas,
           data,
           model_call,
           ndraws = 1000,
           save.files = TRUE,
           path = ".",
           suffix = NULL,
           cores = 1,
           name = NULL) {
    # if (dir.exists(file.path(path, name))){
    #   cat("Directory already existed. Attempting to fit missing models\n")
    # cat("Fitting models (step 1 out of 2) ...")
    # } else
    # dir.create(path = file.path(path, name))
    
    cat("Fitting models (step 1 out of 2) ...")
    fit_list <-
      pblapply_brmsish_int(X = formulas, cl = cores, function(y) {
        
        # make file name without special characters
        mod_name <-
          paste0(remove_special_chars(as.character(y)), ".RDS")
        
        if (save.files &
            !file.exists(file.path(path, mod_name))) {
          cat("Fitting", y, "\n")
          mc <-
            gsub(pattern = "formula",
                 replacement = as.character(y),
                 x = model_call)
          
          mc <- parse(text = mc)
          
          fit <- eval(mc)
          
          if (save.files)
            saveRDS(fit, file = file.path(path, mod_name))
          
        } else {
          cat("Reading", y, "(already existed)\n")
          fit <- readRDS(file.path(path, mod_name))
        }
        return(fit)
      })
    
    if (length(formulas) > 1){
    cat("Averaging models (step 2 out of 2) ...")
    average_call <-
      parse(text = paste("pa(",
                         paste(
                           paste0("fit_list[[", seq_along(fit_list), "]]"), collapse = ", "
                         ),
                         ", ndraws = ",
                         ndraws,
                         ")"))
    
    # Evaluate the expression to create the call object
    average_eval <- eval(average_call)
    
    # add formula as attribute
    attributes(average_eval)$averaged_fit_formulas <- formulas
    
    rds_name <- if (is.null(suffix)) file.path(path, paste0(name, ".RDS")) else
      file.path(path, paste0(suffix, "_", name, ".RDS"))
    
    if (save.files)
      saveRDS(average_eval, file = rds_name)
    
    # return draws from average models
    return(average_eval)
    } else 
    cat("No model averaging conducted as a single formula was supplied")
  }

to_change_percentage <- function(x){
  
  (exp(x) - 1) * 100
  
}

1 Prepare data

1.1 Read data

Code 
clim_dat_2020 <- read_excel("./data/datos_metereologicos/Estacion_met/Dat_met_estacion_2022-01-11.xlsx",sheet = "2020")

clim_dat_2019 <- read_excel("./data/datos_metereologicos/Estacion_met/Dat_met_estacion_2022-01-11.xlsx",sheet = "2019")

clim_dat <- rbind(clim_dat_2019, clim_dat_2020)

clim_dat <- clim_dat[, c("filename", "Año", "Mes", "Día", "Hora", "Temp (°C)", "Humedad Relat.", "Precipitación")]

names(clim_dat) <- c("filename", "year", "month", "day", "hour", "temp", "HR", "rain")

clim_dat <- aggregate(cbind(rain, temp, HR) ~ filename + year + month + day + hour, clim_dat, mean)

clim_dat$year <- clim_dat$year + 2000

clim_dat$date <- as.Date(paste(clim_dat$year, clim_dat$month, clim_dat$day, sep = "-"))


clim_dat$date_hour <- paste(gsub("-", "", clim_dat$date), clim_dat$hour, sep = "-")

call_dat <- read.csv("./data/processed/call_rate_per_date_time_and_site.csv")

# remove 5 pm data and keep only form Sukia
call_dat <- call_dat[call_dat$hour != 17, ]
call_dat <- call_dat[call_dat$site != "LAGCHIMU", ]


call_dat$date_hour <- paste(sapply(as.character(call_dat$site_date_hour), function(x) strsplit(x, "_")[[1]][2]), call_dat$hour, sep = "-") 


call_dat$temp <- sapply(1:nrow(call_dat), function(x){
  y <- clim_dat$temp[clim_dat$date_hour == call_dat$date_hour[x]]
  
  if (length(y) < 1) y <- NA
  
  return(y)
  }
  )


call_dat$HR <- sapply(1:nrow(call_dat), function(x){
  y <- clim_dat$HR[clim_dat$date_hour == call_dat$date_hour[x]]
  
  if (length(y) < 1) y <- NA
  
  return(y)
  }
  )


call_dat$rain <- sapply(1:nrow(call_dat), function(x){
  y <- clim_dat$rain[clim_dat$date_hour == call_dat$date_hour[x]]
  
  if (length(y) < 1) y <- NA
  
  return(y)
  }
  )

# proportion of acoustic data with climatic data
sum(call_dat$date_hour %in% clim_dat$date_hour) / nrow(call_dat)
sum(!is.na(call_dat$temp)) / nrow(call_dat)


call_dat <- call_dat[!is.na(call_dat$temp), ]

call_dat$day <- as.numeric(substr(sapply(as.character(call_dat$site_date_hour), function(x) strsplit(x, "_")[[1]][2]), 7, 8))

call_dat$date <- as.Date(paste(call_dat$year, call_dat$month, call_dat$day, sep = "-"))

call_dat$moon.date <- ifelse(call_dat$hour < 12, as.Date(call_dat$date - 1), as.Date(call_dat$date))

call_dat$moon.date <- as.Date(call_dat$moon.date, origin = "1970-01-02")

## add moon
call_dat$moonlight <- lunar.illumination(call_dat$moon.date, shift = -6)


call_dat$date_hour_min <- strptime(paste(paste(call_dat$year, call_dat$month, call_dat$day, sep = "-"), paste(call_dat$hour, "00", sep = ":")), format="%Y-%m-%d  %H:%M")

call_dat$hour_diff <- as.numeric(call_dat$date_hour_min - min(call_dat$date_hour_min)) / 3600


call_dat$rain_24 <- sapply(1:nrow(call_dat), function(x) sum(clim_dat$rain[strptime(clim_dat$date, format="%Y-%m-%d") == (strptime(call_dat$date[x], format="%Y-%m-%d") - 60 * 60 * 24)]))

call_dat$rain_48 <- sapply(1:nrow(call_dat), function(x) sum(clim_dat$rain[strptime(clim_dat$date, format="%Y-%m-%d") == (strptime(call_dat$date[x], format="%Y-%m-%d") - 60 * 60 * 48)]))

clim_dat$date_hour_min <- strptime(paste(paste(clim_dat$year, clim_dat$month, clim_dat$day, sep = "-"), paste(clim_dat$hour, "00", sep = ":")), format="%Y-%m-%d %H:%M")

clim_dat$hour_diff <- as.numeric(clim_dat$date_hour_min - min(call_dat$date_hour_min)) / 3600

# call_dat$prev_temp <- sapply(1:nrow(call_dat), function(x){
#   
#   # if(call_dat$hour_diff[x] < 48) 
#   #   pt <- NA else
#       pt <- mean(clim_dat$temp[clim_dat$hour_diff %in% (call_dat$hour_diff[x] - 48):(call_dat$hour_diff[x] - 24)])
#     
#     return(pt)
# })

write.csv(call_dat, "./data/processed/acoustic_and_climatic_data_by_hour.csv")

2 Descriptive stats

Code 
call_dat_site <- read.csv("./data/processed/call_rate_per_date_time_and_site.csv")
  • Total number of recordings: 9681
  • Total recordings per site:
Code 
agg_recs <- aggregate(rec_time ~ site, data = call_dat_site, length)
names(agg_recs)[1:2] <- c("site", "rec_count")

# print table as kable  
kb <-kable(agg_recs, row.names = TRUE, digits = 3) 
  
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
  
print(kb)
site rec_count
1 LAGCHIMU 5127
2 SUKIA 4554

  • Total recording time: 3224

  • Total recording time per site:

Code 
agg_recs <- aggregate(rec_time ~ site, data = call_dat_site, sum)
names(agg_recs)[1:2] <- c("site", "recording_time")

agg_recs$recording_time <- round((agg_recs$recording_time) / 60)

# print table as kable  
kb <-kable(agg_recs, row.names = TRUE, digits = 3) 
  
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
  
print(kb)
site recording_time
1 LAGCHIMU 1707
2 SUKIA 1517

  • Total detections: 540203

  • Total detections per site:

Code 
agg_recs <- aggregate(n_call ~ site, data = call_dat_site, sum)
names(agg_recs)[1:2] <- c("calls", "count")

# print table as kable  
kb <-kable(agg_recs, row.names = TRUE, digits = 3) 
  
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
  
print(kb)
calls count
1 LAGCHIMU 169485
2 SUKIA 370718

  • Call rate: 18.598587

  • Call rate per site:

Code 
agg_recs <- aggregate(call_rate ~ site, data = call_dat_site, mean)
agg_recs$sd <- aggregate(call_rate ~ site, data = call_dat_site, FUN = stats::sd)[,2]

names(agg_recs)[1:3] <- c("site", "call_rate", "sd")

print_data_frame(agg_recs)
site call_rate sd
1 LAGCHIMU 11.018 18.039
2 SUKIA 27.133 87.778
Code 
call_rate_hour <- read.csv("./data/processed/acoustic_and_climatic_data_by_hour.csv")


agg <- aggregate(cbind(temp, prev_temp, HR, rain, rain_24, rain_48, moonlight)~1, call_rate_hour, function(x) round(c(mean(x), sd(x), min(x), max(x)), 3))

agg <- as.data.frame(matrix(unlist(agg), ncol = 4, byrow = TRUE, dimnames = list(c("Temperature", "Previous temperature", "Relative humidity", "Night rain", "Rain 24 hours", "Rain 48 hours", "Moonlight"),c("mean", "sd", "min", "max"))))
# print table as kable  
kb <-kable(agg, row.names = TRUE, digits = 3) 
  
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
  
print(kb)
mean sd min max
Temperature 24.229 2.126 19.472 31.763
Previous temperature 24.097 0.988 20.005 27.111
Relative humidity 90.056 8.464 55.158 99.940
Night rain 0.079 0.457 0.000 10.417
Rain 24 hours 1.441 3.161 0.000 25.620
Rain 48 hours 1.453 3.170 0.000 25.620
Moonlight 0.495 0.355 0.000 1.000

Mean and SD temperature: 24.229179 (2.126374)

Mean previous temperature: 24.097073

Mean temperature: 24.229179

Mean cumulative rain per hour: 0.079232

Mean cumulative rain per hour previous 24 hours: 1.44064

Mean daily cumulative rain per hour previous 48 hours: 1.45317

3 Directed acyclical graphs (DAGs)

Code 
coords <- list(
  x = c(
    sc_rain = -0.4,
    evotranspiration = 0.5,
    sc_prev_rain = 0.7,
    sc_temp = -0.8,
    sc_HR = 0,
    n_call = 0,
    sc_moonlight = 0.3,
    hour = -0.5
  ),
  y = c(
    sc_rain = 0.4,
    evotranspiration = 0.3,
    sc_prev_rain = -0.5,
    sc_temp = 0,
    climate = 1,
    sc_HR = -0.6,
    n_call = 0,
    sc_moonlight = 1,
    hour = 0.9
  )
)

# sc_temp + sc_HR + sc_moonlight + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour
# sc_temp = temp y meanT = prev_temp

dag_l <- dagify(sc_rain ~ evotranspiration,
                sc_prev_rain ~ evotranspiration, 
               sc_temp ~ climate,
               sc_temp ~ sc_rain,
               sc_HR ~ sc_rain,
               n_call ~ sc_HR,
               n_call ~ hour,
               n_call ~ sc_moonlight,
               sc_moonlight ~ hour,
               sc_temp ~ hour,
               sc_HR ~ sc_temp,
               sc_HR ~ sc_prev_rain,
               sc_HR ~ sc_rain,
               n_call ~ sc_temp,
               n_call ~ sc_prev_rain,
               n_call ~  sc_rain,
               labels = c("n_call" = "Call\nrate", "sc_HR" = "Relative\nhumidity","sc_rain" = "Current\nrain", "sc_prev_rain" = "Prior\nrain", "sc_moonlight" = "Moonlight", "hour" = "Earth\nrotation", "sc_temp" = "Tempera-\nture", "evotranspiration" = "Evotrans-\npiration", "climate" = "Climate", latent = c("evotranspiration", "climate"), outcome = "n_call"), coords = coords)

tidy_dag <- tidy_dagitty(dag_l)
tidy_dag$data$type <- ifelse(is.na(tidy_dag$data$to), "outcome", "predictor")
tidy_dag$data$type[tidy_dag$data$name %in% c("evotranspiration", "climate")] <- "latent" 


dat <- tidy_dag$data
shorten_distance <- c(0.07, 0.07)
dat$slope <- (dat$yend - dat$y) / (dat$xend - dat$x)
distance <- sqrt((dat$xend - dat$x)^2 + (dat$yend - dat$y)^2)
proportion <- shorten_distance[1]/distance
dat$xend <- (1 - proportion/2) * dat$xend + (proportion/2 * dat$x)
dat$yend <- (1 - proportion/2) * dat$yend + (proportion/2*dat$y)
proportion <- shorten_distance[2]/distance
dat$xstart <- (1-proportion/2)*(dat$x - dat$xend) + dat$xend
dat$ystart <- (1-proportion/2)*(dat$y - dat$yend) + dat$yend

tidy_dag$data <- dat


basic_dag <- ggplot(tidy_dag, aes(
  x = x,
  y = y,
  xend = xend,
  yend = yend
)) +
  scale_color_viridis_d(begin = 0.2, end = 0.8, alpha = 0.5) +
  geom_dag_text(color = "black", aes(label = label, color = label)) + labs(color = "Type") +
  theme_dag() + theme(legend.position = "bottom") + guides(colour = guide_legend(override.aes = list(size =
                                                                                                       10))) +  geom_dag_point(aes(color = type), size = 30) + expand_limits(y = c(-0.67, 1.1))  +
  geom_dag_edges_fan(
    edge_color = viridis(10, alpha = 0.4)[2],
    arrow = grid::arrow(length = grid::unit(10, "pt"), type = "closed"),
    aes(
      x = xstart,
      y = ystart,
      xend = xend,
      yend = yend
    )
  )


ggsave(plot = basic_dag, filename = "./output/dag_300dpi_v02.png", width = 9, height = 5, dpi = 300, bg = "white")


dag_24 <- dagify(sc_rain ~ evotranspiration,
                sc_rain_24 ~ evotranspiration, 
                sc_temp ~ climate,
               sc_temp ~ sc_rain,   
               sc_HR ~ sc_rain,
               n_call ~ sc_HR,
               n_call ~ hour,
               n_call ~ sc_moonlight,
               sc_moonlight ~ hour,
               sc_temp ~ hour,
               sc_HR ~ sc_temp,
               sc_HR ~ sc_rain_24,
               sc_HR ~ sc_rain,
               n_call ~ sc_temp,
               n_call ~ sc_rain_24,
               n_call ~  sc_rain,
               labels = c("n_call" = "Call rate", "sc_HR" = "Relative humidity","sc_rain" = "Night Rain", "sc_rain_24" = "Previous Rain", "sc_moonlight" = "Moonlight", "hour" = "Earth rotation", "sc_temp" = "Temperature", "evotranspiration" = "Evotranspiration", "climate" = "Climate", latent = c("evotranspiration", "climate"), outcome = "n_call"))

dag_48 <- dagify(sc_rain ~ evotranspiration,
                sc_rain_48 ~ evotranspiration, 
               sc_temp ~ climate,
               sc_temp ~ sc_rain,
               sc_HR ~ sc_rain,
               n_call ~ sc_HR,
               n_call ~ hour,
               n_call ~ sc_moonlight,
               sc_moonlight ~ hour,
               sc_temp ~ hour,
               sc_HR ~ sc_temp,
               sc_HR ~ sc_rain_48,
               sc_HR ~ sc_rain,
               n_call ~ sc_temp,
               n_call ~ sc_rain_48,
               n_call ~  sc_rain,
               labels = c("n_call" = "Call rate", "sc_HR" = "Relative humidity","sc_rain" = "Night Rain", "sc_rain_48" = "Previous Rain", "sc_moonlight" = "Moonlight", "hour" = "Earth rotation", "sc_temp" = "Temperature", "evotranspiration" = "Evotranspiration", "climate" = "Climate", latent = c("evotranspiration", "climate"), outcome = "n_call"))

4 Bayesian regression models

4.1 Scale variables and set model parameters

Code 
call_rate_hour <- read.csv("./data/processed/acoustic_and_climatic_data_by_hour.csv")

# make hour a factor
call_rate_hour$hour <- factor(call_rate_hour$hour)

# scale and mean-center
call_rate_hour$sc_temp <- scale(call_rate_hour$temp)
call_rate_hour$sc_HR <- scale(call_rate_hour$HR)
call_rate_hour$sc_rain <- scale(call_rate_hour$rain)
call_rate_hour$sc_rain_24 <- scale(call_rate_hour$rain_24)
call_rate_hour$sc_rain_48 <- scale(call_rate_hour$rain_48)
call_rate_hour$sc_moonlight <- scale(call_rate_hour$moonlight)

priors <- c(prior(normal(0, 4), class = "b"))
chains <- 4
iter <- 10000

4.2 Fit all models

Code 
param_grid <- expand.grid(
  dag = c("dag_24", "dag_48"), 
  exposure = c("sc_temp", "sc_HR", "sc_moonlight", "sc_rain", "sc_rain_24", "sc_rain_48"),
  effect = c("total", "direct"), stringsAsFactors = FALSE
  )

param_grid$name <- apply(param_grid, 1, paste, collapse = "-")

# remove wrong dags for previous rain
param_grid <- param_grid[!(param_grid$dag == "dag_24" & param_grid$exposure == "sc_rain_48") & !(param_grid$dag == "dag_48" & param_grid$exposure == "sc_rain_24"), ]


adjustment_sets_list <- pblapply(seq_len(nrow(param_grid)), cl = 1, function(x){
  
forms <- adjustment_set_formulas(
  dag = if (param_grid$dag[x] == "dag_24") dag_24 else dag_48,
  type = if (param_grid$effect[x] == "total") "all" else "minimal",
  exposure = param_grid$exposure[x],
  outcome = "n_call",
  effect = param_grid$effect[x],
  required_variable = "hour",
  formula_parts = c(
    "n_call | resp_rate(rec_time)",
    "+ ar(p = 2, time = hour_diff, gr = hour)"
  ),
  latent = c("evotranspiration", "climate"),
  remove = "hour",
  plot = FALSE
)

return(forms)
})



names(adjustment_sets_list) <- param_grid$name

param_grid$model <-sapply(seq_len(nrow(param_grid)), function(x) {
  if (param_grid$effect[x] == "direct")
  adjustment_sets_list[[which(names(adjustment_sets_list) == param_grid$name[x])]] else
    NA
  }) 

param_grid$model <- unlist(
param_grid$model)
param_grid <- as.data.frame(param_grid)


param_grid$model[!is.na(param_grid$model)] <- remove_special_chars(param_grid$model[!is.na(param_grid$model)] ) 
param_grid$model <- c("total_effect_temperature_with_rain_24", "total_effect_temperature_with_rain_48",  "total_effect_humidity_with_rain_24",  "total_effect_humidity_with_rain_48", "total_effect_moon_with_rain_24", "total_effect_moon_with_rain_48", "total_effect_rain_with_rain_24", "total_effect_rain_with_rain_48", "total_effect_previous_rain_24", "total_effect_previous_rain_48", param_grid$model[!is.na(param_grid$model)])

param_grid$exposure.name <- param_grid$exposure
param_grid$exposure.name[grep("temp", param_grid$exposure.name)] <- "Temperature"
param_grid$exposure.name[grep("HR", param_grid$exposure.name)] <- "Relative humidity"
param_grid$exposure.name[grep("moon", param_grid$exposure.name)] <- "Moonlight"
param_grid$exposure.name[grep("rain$", param_grid$exposure.name)] <- "Current rain"
param_grid$exposure.name[grep("rain_24", param_grid$exposure.name)] <- "Previous rain (24h)"
param_grid$exposure.name[grep("rain_48", param_grid$exposure.name)] <- "Previous rain (48h)"

table(param_grid$exposure.name)

write.csv(x = param_grid, file = "./data/processed/direct_and_total_effect_model_data_frame.csv", row.names = FALSE)

direct_adjustment_sets_list <- adjustment_sets_list[grep("direct", names(adjustment_sets_list))]


for(i in seq_along(direct_adjustment_sets_list))
pa_comb_mod <-
averaged_model(
    formulas = direct_adjustment_sets_list[[i]],
    data = call_rate_hour,
    suffix = "direct",
    model_call = "brm(formula, data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(), prior = priors, backend = 'cmdstanr')",
    save.files = TRUE,
    path = "./data/processed/averaged_models", 
    # name = "temperature_with_rain_24",
    cores = 1 
  )


model_call = "brm(formula, data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(), prior = priors, backend = 'cmdstanr')"
formulas <- unlist(direct_adjustment_sets_list)
path <-  "./data/processed/single_models"

fit_list <- pblapply_brmsish_int(X = formulas, cl = 1, function(y) {
        # make file name without special characters
        mod_name <-
          paste0(remove_special_chars(as.character(y)), ".RDS")
        
        if (!file.exists(file.path(path, mod_name))) {
          cat("Fitting", y, "\n")
          mc <-
            gsub(pattern = "formula",
                 replacement = as.character(y),
                 x = model_call)
          
          mc <- parse(text = mc)
          
          fit <- eval(mc)
          
          if (save.files)
            saveRDS(fit, file = file.path(path, mod_name))
          
        } 
      })

5 Results

Code 
param_grid <- read.csv(file = "./data/processed/direct_and_total_effect_model_data_frame.csv")

param_grid$files <- file.path("./data/processed/averaged_models", paste0(param_grid$model, ".RDS"))


for(i in unique(param_grid$exposure.name)){
  Y <- param_grid[param_grid$exposure.name == i, ]
  
  cat(paste("\n##", i), "\n")
  cat("\n### Direct effects\n")
  for(e in which(Y$effect == "direct")){
    if (grepl("24", Y$model[e]))
      cat("\n#### 24 hour previous rain model:\n") else
              cat("\n#### 48 hour previous rain model:\n")
    extended_summary(
    read.file = Y$files[e],
    highlight = TRUE,
    remove.intercepts = TRUE, 
    print.name = FALSE
    )
    cat("\n")
    }

 cat("\n### Total effect\n")
  for(w in which(Y$effect == "total")){
       if (grepl("24", Y$files[w]))
      cat("\n#### 24 hour previous rain model:\n") else
              cat("\n#### 48 hour previous rain model:\n")
    
    draws <- readRDS(Y$files[w])

draw_extended_summary(
  draws,
  highlight = TRUE,
  remove.intercepts = TRUE
)
  cat("\n")
 cat("\n##### Summary of single models:\n")
    
  # print summary
  print(readRDS(gsub("\\.RDS", "_fit_summary.RDS", Y$files[w])))
  }
    cat("\n")
  }

5.1 Temperature

5.1.1 Direct effects

5.1.1.1 24 hour previous rain model:

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 13808 14519.5 1687235322
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_temp 0.577 0.498 0.656 1 13808.0 14685.6
b_sc_HR 0.314 0.242 0.387 1 18992.0 14519.5
b_sc_rain 0.037 0.002 0.074 1 23964.4 15266.3
b_sc_rain_24 0.096 0.057 0.135 1 23319.6 15914.6

5.1.1.2 48 hour previous rain model:

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 9460.49 12432.4 293904519
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_temp 0.565 0.487 0.643 1 9460.49 12432.4
b_sc_HR 0.307 0.235 0.379 1 11994.69 13312.6
b_sc_rain 0.029 -0.007 0.066 1 17566.34 16050.0
b_sc_rain_48 0.006 -0.032 0.044 1 17289.44 16222.0

5.1.2 Total effect

5.1.2.1 24 hour previous rain model:

Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_temp 0.308 0.258 0.361 1.043 44.190 893.839
b_sc_rain 0.041 0.004 0.078 1.026 152.336 830.415

5.1.2.1.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2542.43 5383.76 629626678
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2238.34 5028.64 2096770826
3 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2320.90 4619.08 1430954384
4 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2292.54 4473.95 1881604582

5.1.2.2 48 hour previous rain model:

Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_temp 0.305 0.258 0.353 1.005 961.292 983.283
b_sc_rain 0.036 -0.002 0.076 1 1113.522 969.564

5.1.2.2.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2529.79 5016.37 1471048663
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2588.68 5844.63 412112948
3 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2622.63 5641.91 41093237
4 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2334.61 4741.38 366114719

5.2 Relative humidity

5.2.1 Direct effects

5.2.1.1 24 hour previous rain model:

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 11039.8 13170.1 1554924087
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_HR 0.314 0.240 0.388 1 13006.8 13679.9
b_sc_rain 0.037 0.002 0.075 1 18235.1 15400.0
b_sc_rain_24 0.096 0.057 0.135 1 16966.8 15476.0
b_sc_temp 0.578 0.499 0.657 1 11039.8 13170.1

5.2.1.2 48 hour previous rain model:

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 8448.99 11500.1 1931567299
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_HR 0.308 0.234 0.380 1 10649.49 12433.6
b_sc_rain 0.029 -0.008 0.067 1 16664.06 15531.5
b_sc_rain_48 0.005 -0.032 0.043 1 15745.61 14835.4
b_sc_temp 0.566 0.487 0.645 1 8448.99 11500.1

5.2.2 Total effect

5.2.2.1 24 hour previous rain model:

Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_HR 0.309 0.240 0.383 1.006 942.064 903.519
b_sc_rain 0.038 0.002 0.076 1.003 995.090 912.254
b_sc_rain_24 0.094 0.057 0.133 1.003 993.441 772.255
b_sc_temp 0.575 0.496 0.655 1.004 1022.397 861.998

5.2.2.1.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_HR + sc_moonlight + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2556.01 5127.35 1982612331
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2182.06 4328.43 1554924087

5.2.2.2 48 hour previous rain model:

Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_HR 0.303 0.230 0.371 1.006 801.527 931.946
b_sc_rain 0.030 -0.008 0.066 1.002 812.238 1071.444
b_sc_rain_48 0.003 -0.034 0.044 1.005 889.913 984.434
b_sc_temp 0.564 0.480 0.644 1.003 821.550 773.876

5.2.2.2.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_HR + sc_moonlight + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2326.55 4696.55 846435536
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2166.28 4627.63 1931567299

5.3 Moonlight

5.3.1 Direct effects

5.3.1.1 48 hour previous rain model:

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 30283.5 15929.8 1624275037
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_moonlight -0.129 -0.193 -0.064 1 30283.5 15929.8

5.3.1.2 48 hour previous rain model:

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 30283.5 15929.8 1624275037
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_moonlight -0.129 -0.193 -0.064 1 30283.5 15929.8

5.3.2 Total effect

5.3.2.1 24 hour previous rain model:

Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_moonlight -0.101 -0.17 -0.039 1.04 40.665 750.448

5.3.2.1.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2463.32 4856.49 1624275037
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2613.63 5435.51 1331671454
3 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2426.99 4881.26 1714115616
4 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2184.13 4557.11 946500998
5 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2386.65 5083.35 2088023938
6 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2571.23 4971.80 1587538237
7 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2319.91 4626.15 1236984196
8 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2197.61 4746.52 1987527846
9 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2353.89 5263.44 924618251
10 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2409.47 4652.12 1063522846
11 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2407.64 4851.64 259660598
12 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2473.17 4731.37 1625659188
13 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2596.10 5292.48 1259975215
14 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 1998.35 4047.25 192228769
15 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2365.76 4743.96 559459570
16 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2432.09 4364.34 1786523665

5.3.2.2 48 hour previous rain model:

Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_moonlight -0.105 -0.166 -0.038 1.02 111.506 659.054

5.3.2.2.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2534.83 5170.40 169376128
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2664.91 5585.82 731739252
3 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2449.93 4958.79 230289957
4 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2283.15 4266.04 1735286727
5 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2486.75 4561.76 1539258099
6 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2525.61 5053.20 1234320719
7 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2481.99 4862.29 183236870
8 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2296.24 4114.83 239469808
9 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2298.87 4495.74 1126617762
10 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2442.38 4594.40 684931243
11 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 1 (5e-05%) 0 2416.00 4930.19 1093486403
12 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2146.30 4395.84 1670597608
13 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2450.19 4442.71 627408872
14 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2513.43 4849.68 1827601939
15 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2101.51 4415.14 104020182
16 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2495.11 5464.93 598459253
17 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2120.85 4487.97 2041361978
18 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2336.51 4816.06 1759350038
19 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 1 (5e-05%) 0 2453.12 5040.07 927628753
20 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2360.45 5250.70 2028572807
21 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2642.43 5077.10 293339431
22 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2342.88 4887.79 1868346384
23 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2515.33 5176.64 552665161
24 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2253.53 4141.73 1764642823
25 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2634.31 5119.83 505424120
26 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2427.59 5073.57 1761124445
27 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2609.57 5134.76 1073467634
28 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2461.43 4814.97 1419285187
29 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2596.36 5583.40 1232772932
30 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2271.50 5053.09 278133969
31 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 1 (5e-05%) 0 2507.89 4315.91 866725413
32 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2126.46 4549.18 1885323046

5.4 Current rain

5.4.1 Direct effects

5.4.1.1 24 hour previous rain model:

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain + sc_HR + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 10228.7 12456.7 1281527469
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain 0.037 0.001 0.075 1 18004.9 15909.5
b_sc_HR 0.313 0.240 0.386 1 13025.4 13472.6
b_sc_rain_24 0.096 0.059 0.134 1 17444.9 15779.0
b_sc_temp 0.577 0.499 0.655 1 10228.7 12456.7

5.4.1.2 48 hour previous rain model:

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain + sc_HR + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 7371.36 11882.6 1129340659
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain 0.029 -0.007 0.066 1 14353.91 14479.2
b_sc_HR 0.307 0.235 0.380 1 8890.14 12279.5
b_sc_rain_48 0.005 -0.032 0.043 1 13843.68 14097.9
b_sc_temp 0.566 0.487 0.644 1 7371.36 11882.6

5.4.2 Total effect

5.4.2.1 24 hour previous rain model:

Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain 0.000 -0.034 0.035 1.001 942.79 942.018
b_sc_rain_24 0.076 0.039 0.113 1 1056.52 982.799

5.4.2.1.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain + sc_moonlight + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2531.75 4835.39 1596946440
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2369.16 4574.59 1687235322

5.4.2.2 48 hour previous rain model:

Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain -0.005 -0.042 0.030 1 1024.455 819.064
b_sc_rain_48 0.011 -0.027 0.049 1.009 932.089 975.233

5.4.2.2.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain + sc_moonlight + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2647.66 5129.91 1485832077
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2675.26 5099.07 167092432

5.5 Previous rain (24h)

5.5.1 Direct effects

5.5.1.1 24 hour previous rain model:

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 11374.2 13445.2 1379661434
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain_24 0.096 0.058 0.135 1 20345.7 15663.3
b_sc_HR 0.314 0.240 0.386 1 14938.4 13445.2
b_sc_rain 0.038 0.001 0.074 1 21502.1 16005.3
b_sc_temp 0.578 0.498 0.656 1 11374.2 13493.9

5.5.2 Total effect

5.5.2.1 24 hour previous rain model:

Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain_24 0.088 0.051 0.127 1.01 883.631 733.445
b_sc_rain 0.043 0.005 0.081 1.001 908.203 446.195

5.5.2.1.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2692.30 5204.02 423216939
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2556.12 4635.90 610868624
3 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 1 (5e-05%) 0 2559.09 4449.50 1988187279
4 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2077.52 4134.53 640249946

5.6 Previous rain (48h)

5.6.1 Direct effects

5.6.1.1 48 hour previous rain model:

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain_48 + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 11112 13056.2 1238572765
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain_48 0.006 -0.031 0.043 1 17211.0 15210.0
b_sc_HR 0.308 0.236 0.380 1 14154.7 13956.7
b_sc_rain 0.029 -0.007 0.066 1 19170.2 15964.2
b_sc_temp 0.566 0.486 0.645 1 11112.0 13056.2

5.6.2 Total effect

5.6.2.1 48 hour previous rain model:

Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain_48 0.005 -0.033 0.043 1.014 711.511 873.701
b_sc_rain 0.036 0.001 0.073 1 857.372 792.928

5.6.2.1.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2447.89 5823.08 1081516016
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 2 (1e-04%) 0 2413.93 4851.86 217409394
3 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2506.36 5173.40 724206553
4 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2641.41 5493.32 509251240
5 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_48 + sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2386.41 4640.41 518390312
6 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_48 + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2575.04 5332.87 76427694

5.7 Combined results with causal inference estimates

Takes the posterior probability distributions from the right causal models

Code 
param_grid <- read.csv(file = "./data/processed/direct_and_total_effect_model_data_frame.csv")
param_grid <- param_grid[param_grid$effect == "direct",]

param_grid$file <- paste0(remove_special_chars(param_grid$model), ".RDS")

rdss_24 <- list.files("./data/processed/averaged_models", pattern = "24.RDS", full.names = TRUE)

combined_draws_list <- lapply(rdss_24, function(x) {
  
  total_draws <- readRDS(x)
  
  exp <-
    attributes((attributes(total_draws)$averaged_fit_formulas))$exposure.formula
  exp <-
    gsub("n_call | resp_rate(rec_time)  ~  ", "", exp, fixed = TRUE)
  exposure <-
    exp <-
    gsub(" + ar(p = 2, time = hour_diff, gr = hour)", "", exp, fixed = TRUE)
  exp <- paste0("b_", exp)
  total_draws <-
    total_draws[, colnames(total_draws) == exp, drop  = FALSE]
  names(total_draws) <- exp
  
  direct_fit_file <-
    param_grid$file[param_grid$exposure == exposure]
  
  direct_fit_file <- direct_fit_file[!duplicated(direct_fit_file)]
  
  if (length(direct_fit_file) > 1)
    direct_fit_file <- grep("24", direct_fit_file, value = TRUE)
  
  direct_fit <-
    readRDS(file = file.path("./data/processed/averaged_models", direct_fit_file))
  direct_draws <-
    posterior::merge_chains(as_draws(direct_fit, variable = exp))
  direct_draws <-
    as.data.frame(thin_draws(direct_draws, thin = length(direct_draws[[1]][[1]])
                             / (nrow(total_draws)))[[1]])
  
  direct_draws$effect <- "direct"
  total_draws$effect <- "total"
  
  draws <- rbind(direct_draws, total_draws)
  
  return(draws)
})

combined_draws <- do.call(cbind, combined_draws_list)
combined_draws <- combined_draws[, c(which(sapply(combined_draws, is.numeric)), ncol(combined_draws))]

combined_draws[,-ncol(combined_draws)] <- to_change_percentage(combined_draws[,-ncol(combined_draws)])


# combined_draws <- as.data.frame(combined_draws)
combined_draws$effect <- ifelse(combined_draws$effect == "direct", "Direct", "Total")


saveRDS(combined_draws, "./data/processed/combined_draws_for_total_and_direct_effects_24h_previous_rain.RDS")
Code 
combined_draws <- readRDS("./data/processed/combined_draws_for_total_and_direct_effects_24h_previous_rain.RDS")

fill_colors <- viridis::mako(10)[c(8, 4)]

gg_dists <- draw_extended_summary(
  draws = combined_draws,
  highlight = TRUE,
  remove.intercepts = TRUE,
  fill = adjustcolor(fill_colors, alpha.f = 0.4),
  by = "effect",
  gsub.pattern = c(
    "b_sc_HR",
    "b_sc_rain$",
    "b_sc_rain_24",
    "b_sc_temp",
    "b_sc_moonlight"
  ),
  gsub.replacement = c(
    "Relative\nhumidity",
    "Current\nrain",
    "Prior\nrain",
    "Temperature",
    "Moonlight"
  ), 
  ylab = "Variable",
  xlab = "Effect size (% of change)"
)
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Direct-Current rain 3.652 -0.002 7.942 1.009 999.221 949.102
Total-Current rain 0.000 -3.373 3.596 1.001 942.790 942.018
Direct-Moonlight -12.018 -17.074 -6.217 1.001 972.135 944.341
Total-Moonlight -9.600 -15.624 -3.860 1.04 40.665 750.448
Direct-Prior rain 10.060 5.793 14.735 1 916.714 983.283
Total-Prior rain 9.189 5.189 13.581 1.01 883.631 733.445
Direct-Relative humidity 37.000 27.013 47.025 1.001 902.884 848.891
Total-Relative humidity 36.187 27.104 46.613 1.006 942.064 903.519
Direct-Temperature 77.591 65.145 92.509 1.002 955.952 875.689
Total-Temperature 36.112 29.497 43.428 1.043 44.190 893.839

Posterior distribution of direct (green) and total (purple) effect sizes of environmental factors on the calling activity of A. lemur. Posterior values were transformed into percentage change to facilitate interpretation. Dots and error bars show the median and 95% uncertainty intervals of the distributions. Solid color distributions correspond to effect sizes in which uncertainty intervals do not include zero. “Prior rain” accounts for the 24 hour period before sampling calling activity.
Code 
gg_dists + ggplot2::scale_fill_manual(values = fill_colors) +
ggplot2::theme(
      axis.ticks.length = ggplot2::unit(0, "pt"),
      plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"),
      legend.position = "inside",
      legend.position.inside = c(0.7, 0.7)
      )

Posterior distribution of direct (green) and total (purple) effect sizes of environmental factors on the calling activity of A. lemur. Posterior values were transformed into percentage change to facilitate interpretation. Dots and error bars show the median and 95% uncertainty intervals of the distributions. Solid color distributions correspond to effect sizes in which uncertainty intervals do not include zero. “Prior rain” accounts for the 24 hour period before sampling calling activity.
Code 
ggsave("./output/posterior_distribution_change_percentage_24h_previous_rain.png", width = 5, height = 3, dpi = 300)
Code 
param_grid <- read.csv(file = "./data/processed/direct_and_total_effect_model_data_frame.csv")
param_grid <- param_grid[param_grid$effect == "direct",]

param_grid$file <- paste0(remove_special_chars(param_grid$model), ".RDS")

rdss_48 <- list.files("./data/processed/averaged_models", pattern = "48.RDS", full.names = TRUE)

combined_draws_list <- lapply(rdss_48, function(x) {

  total_draws <- readRDS(x)
  
  exp <-
    attributes((attributes(total_draws)$averaged_fit_formulas))$exposure.formula
  exp <-
    gsub("n_call | resp_rate(rec_time)  ~  ", "", exp, fixed = TRUE)
  exposure <-
    exp <-
    gsub(" + ar(p = 2, time = hour_diff, gr = hour)", "", exp, fixed = TRUE)
  exp <- paste0("b_", exp)
  total_draws <-
    total_draws[, colnames(total_draws) == exp, drop  = FALSE]
  names(total_draws) <- exp
  
  direct_fit_file <-
    param_grid$file[param_grid$exposure == exposure]
  
  direct_fit_file <- direct_fit_file[!duplicated(direct_fit_file)]
  
  if (length(direct_fit_file) > 1)
    direct_fit_file <- grep("48", direct_fit_file, value = TRUE)
  
  direct_fit <-
    readRDS(file = file.path("./data/processed/averaged_models", direct_fit_file))
  direct_draws <-
    posterior::merge_chains(as_draws(direct_fit, variable = exp))
  direct_draws <-
    as.data.frame(thin_draws(direct_draws, thin = length(direct_draws[[1]][[1]])
                             / (nrow(total_draws)))[[1]])
  
  direct_draws$effect <- "direct"
  total_draws$effect <- "total"
  
  draws <- rbind(direct_draws, total_draws)
  
  return(draws)
})

combined_draws <- do.call(cbind, combined_draws_list)
combined_draws <- combined_draws[, c(which(sapply(combined_draws, is.numeric)), ncol(combined_draws))]

combined_draws[,-ncol(combined_draws)] <- to_change_percentage(combined_draws[,-ncol(combined_draws)])


# combined_draws <- as.data.frame(combined_draws)
combined_draws$effect <- ifelse(combined_draws$effect == "direct", "Direct", "Total") 

saveRDS(combined_draws, "./data/processed/combined_draws_for_total_and_direct_effects_48h_previous_rain.RDS")
Code 
combined_draws <- readRDS("./data/processed/combined_draws_for_total_and_direct_effects_48h_previous_rain.RDS")


gg_dists <- draw_extended_summary(
  draws = combined_draws,
  highlight = TRUE,
  remove.intercepts = TRUE,
  fill = adjustcolor(fill_colors, alpha.f = 0.4),
  by = "effect",
  gsub.pattern = c(
    "b_sc_HR",
    "b_sc_rain$",
    "b_sc_rain_48",
    "b_sc_temp",
    "b_sc_moonlight"
  ),
  gsub.replacement = c(
    "Relative\nhumidity",
    "Current\nrain",
    "Prior\nrain",
    "Temperature",
    "Moonlight"
  ), 
  ylab = "Variable",
  xlab = "Effect size (% of change)"
)
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Direct-Current rain 2.940 -0.477 6.832 1 933.292 847.036
Total-Current rain -0.544 -4.093 3.071 1 1024.455 819.064
Direct-Moonlight -12.018 -17.074 -6.217 1.001 972.135 944.341
Total-Moonlight -9.925 -15.276 -3.765 1.02 111.506 659.054
Direct-Prior rain 0.623 -2.979 4.124 1.003 857.680 1009.747
Total-Prior rain 0.540 -3.268 4.427 1.014 711.511 873.701
Direct-Relative humidity 36.477 26.300 46.337 0.999 1043.482 878.337
Total-Relative humidity 35.405 25.805 44.967 1.006 801.527 931.946
Direct-Temperature 76.199 63.368 89.305 1.004 1009.050 835.503
Total-Temperature 35.603 29.387 42.312 1.005 961.292 983.283

Code 
gg_dists + ggplot2::scale_fill_manual(values = fill_colors) +
ggplot2::theme(
      axis.ticks.length = ggplot2::unit(0, "pt"),
      plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"),
      legend.position = "inside",
      legend.position.inside = c(0.7, 0.7)
      )

Code 
ggsave("./output/posterior_distribution_change_percentage_48h_previous_rain.png", width = 5, height = 3, dpi = 300)

Takeaways

  • Variation in call activity strongly linked to environmental variation

Session information

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] pbapply_1.7-2      tidybayes_3.0.6    ggdag_0.2.12       dagitty_0.3-4     
 [5] ohun_1.0.1         warbleR_1.1.30     NatureSounds_1.0.4 seewave_2.2.3     
 [9] tuneR_1.4.7        brmsish_1.0.0      lunar_0.2-1        ggplot2_3.5.1     
[13] ggdist_3.3.2       knitr_1.47         kableExtra_1.4.0   HDInterval_0.2.4  
[17] readxl_1.4.3       posterior_1.5.0    cowplot_1.1.3      brms_2.21.0       
[21] Rcpp_1.0.12        viridis_0.6.5      viridisLite_0.4.2  remotes_2.5.0     

loaded via a namespace (and not attached):
  [1] tensorA_0.36.2.1     rstudioapi_0.16.0    jsonlite_1.8.8      
  [4] magrittr_2.0.3       TH.data_1.1-2        sketchy_1.0.3       
  [7] estimability_1.5.1   farver_2.1.2         rmarkdown_2.27      
 [10] ragg_1.3.2           vctrs_0.6.5          memoise_2.0.1       
 [13] RCurl_1.98-1.14      htmltools_0.5.8.1    distributional_0.4.0
 [16] curl_5.2.1           signal_1.8-1         cellranger_1.1.0    
 [19] StanHeaders_2.32.9   KernSmooth_2.23-24   htmlwidgets_1.6.4   
 [22] plyr_1.8.9           sandwich_3.1-0       testthat_3.2.1.1    
 [25] cachem_1.1.0         emmeans_1.10.2       zoo_1.8-12          
 [28] igraph_2.0.3         lifecycle_1.0.4      pkgconfig_2.0.3     
 [31] Matrix_1.7-0         R6_2.5.1             fastmap_1.2.0       
 [34] digest_0.6.36        colorspace_2.1-0     textshaping_0.4.0   
 [37] labeling_0.4.3       fansi_1.0.6          polyclip_1.10-6     
 [40] abind_1.4-5          compiler_4.4.1       proxy_0.4-27        
 [43] withr_3.0.0          backports_1.5.0      inline_0.3.19       
 [46] DBI_1.2.3            QuickJSR_1.2.2       pkgbuild_1.4.4      
 [49] highr_0.11           ggforce_0.4.2        MASS_7.3-61         
 [52] rjson_0.2.21         classInt_0.4-10      loo_2.7.0           
 [55] tools_4.4.1          units_0.8-5          ape_5.8             
 [58] fftw_1.0-8           glue_1.7.0           nlme_3.1-165        
 [61] grid_4.4.1           sf_1.0-6             checkmate_2.3.1     
 [64] reshape2_1.4.4       generics_0.1.3       diffobj_0.3.5       
 [67] gtable_0.3.5         class_7.3-22         tidyr_1.3.1         
 [70] tidygraph_1.3.1      xml2_1.3.6           utf8_1.2.4          
 [73] ggrepel_0.9.5        pillar_1.9.0         stringr_1.5.1       
 [76] splines_4.4.1        dplyr_1.1.4          tweenr_2.0.3        
 [79] lattice_0.22-6       survival_3.7-0       tidyselect_1.2.1    
 [82] arrayhelpers_1.1-0   gridExtra_2.3        V8_4.4.2            
 [85] svglite_2.1.3        stats4_4.4.1         xfun_0.45           
 [88] graphlayouts_1.1.1   bridgesampling_1.1-2 brio_1.1.5          
 [91] matrixStats_1.3.0    rstan_2.32.6         stringi_1.8.4       
 [94] yaml_2.3.8           boot_1.3-30          xaringanExtra_0.8.0 
 [97] evaluate_0.24.0      codetools_0.2-20     dtw_1.23-1          
[100] ggraph_2.2.1         tibble_3.2.1         cli_3.6.3           
[103] RcppParallel_5.1.7   xtable_1.8-4         systemfonts_1.1.0   
[106] munsell_0.5.1        coda_0.19-4.1        svUnit_1.0.6        
[109] parallel_4.4.1       rstantools_2.4.0     bayesplot_1.11.1    
[112] Brobdingnag_1.2-9    bitops_1.0-7         mvtnorm_1.2-5       
[115] scales_1.3.0         e1071_1.7-14         packrat_0.9.2       
[118] purrr_1.0.2          crayon_1.5.3         rlang_1.1.4         
[121] multcomp_1.4-25